Marker identification and processing in x-ray images

ABSTRACT

A mechanism for rapidly detecting and localizing external markers placed on a patient in projection images. Embodiments of the invention allow the markers to be detected even in the presence of dense surrounding anatomy and extensive patient motion. Once the positions of the marker points on the projection images are extracted, the marker points can be used to perform marker-based patient motion detection. Embodiments of the invention can also be used outside of motion correction, such as for scanner calibration, automatic cephalometric measurements, and quality control assessment.

BACKGROUND

Embodiments of the invention relate to image processing. In particular,embodiments provide methods and systems for detecting, localizing, andclassifying marker-associated points located in x-ray images.

SUMMARY

Modern imaging systems acquire many images during a single procedure orscan. The images are combined using a process known as reconstruction toproduce one or more images representing finite-width “slices” through apatient. Because the resulting images are often analyzed in preparationfor medical or dental procedures or to properly calibrate the imagingequipment, it is important that the images be as clear as possible.However, patient movement during the scan can cause errors or blurringin the image.

Markers are often placed on a patient before an imaging procedure. Themarkers appear on the resulting image as marker points and are used toidentify particular features in the image. The marker points can also beused to correct the images for patient movement. In a typical motioncorrection scenario, a relatively small number of markers (e.g., seven)must be identified in each of the hundreds (e.g., 600) or more imagescaptured during the scan. Thus, using markers to detect and correctpatient motion is practical only if the thousands of marker points(e.g., 4,200 marker points using the example above) can be automaticallyidentified in the images. Marker point identification is furthercomplicated by the non-linear nature of x-ray imaging and the relativelylarge dynamic range of x-ray images.

Yet another complication is that marker point identification and theresulting motion correction needs to be performed quickly because itmust be completed before reconstruction of the images can begin.Existing methods for identifying marker-like points on images rely onthe general circular nature of these points and employ pointidentification techniques such as image gradient processing or somevariation of a Hough transform, which can take as long as 30 seconds perimage to complete. Even the simplest methods can take 1/10 of a secondor more per image and often take as long as five seconds per image.

The accuracy of these motion correction methods is directly related tothe accuracy in extracting marker point locations. For example, in mostcorrection methods the markers must be localized to a precision of ¼ to⅛ of a pixel. This degree of accuracy must be accomplished in thepresence of a wide dynamic range of background values, a non-lineardependence on image background, and noise. Many existing correctionmethods do not properly deal with the non-linear nature of the imagebackground, which leads to positional errors at image edges. Theseexisting methods also do not work well when markers in images arerelatively close to each other.

Once marker points are identified in an image they also need to bemapped back to the physical markers from which they originated. Thephysical markers are uniquely identified by their three-dimensional(“3D”) positional coordinates. Thus, to make lists of points from eachprojection image associated with particular markers, one step is toidentify each marker's physical location. This can be done byidentifying the 3D location for each marker that is most consistent withhow the position of the points produced by the marker changes fromprojection image to projection image as the gantry rotates. However,each marker's associated point locations may vary significantly fromtheir expected locations because of patient movement or unplanned systemor equipment movement. Therefore, the accuracy of identifying a marker's3-dimensional coordinates suffers in the presence of noise or patientmotion or when multiple markers are positioned close to each other.Existing techniques for identifying a marker's 3D location may use aposition-dependent electro-magnetic field that is detected by a remotesensor and is converted to a position. In addition to the complexity ofthis transmitter and sensor system, the technology used by the systemdepends on multiple sensors to determine a 3-dimensional position of amarker. Other localization techniques involve optical technologies thatdepend on knowledge of a fixed relationship between multiple detectedsites, which increases the complexity of the technology.

Even after the physical markers associated with the image points areidentified and a position is assigned to each, the problem still existsof assigning each point on each image to its source marker. Before thisis done, all the points in all the images exist only as a single list.Point assignment (marker tracking) attempts to associate each markerpoint on each image with the marker that actually generated it. Afterpoint assignment is complete, each marker will possess a separate listof the marker points associated with that marker. Each marker's separatelist may have only one or zero marker points at each projection image.It is important that the correct correspondence between markers andtheir associated marker points be maintained, even in the presence ofmotion and measurement error.

Embodiments provide systems and methods for rapidly identifyingmarker-associated points in an image, which can then be used in motioncorrection systems. One method includes processing each projection imageto reduce the intensity of all objects but circular objects (i.e., theshape of the marker points), identifying local maxima in each image,growing the local maxima to create candidate regions, and applyingcircularity and projection continuity criteria to each candidate region.The method uses a sequence of applied simple mean filters, differences,and division/subtraction operations to reduce the intensity of all theobjects but circular objects, in a manner which reduces the total imageprocessing time by at least an order of a magnitude compared to othermethods. One method includes obtaining an image, applying a 2-D highpass filter to create a background suppressed image, and extractingmarker associated points from the background image.

Embodiments also provide systems and methods for performing accuratesub-pixel localization of marker point locations in an image. One methodincludes independently evaluating the position of each marker pointidentified in an image in each of the two dimensions. This evaluation isdone by determining a background profile of each marker point bysampling one or both sides of a marker point. Taking a ratio of abaseline profile and the marker point profile corrects for backgroundeffects while addressing the non-linear nature of x-ray attenuation.Residual background offset is removed by subtracting the average valueof tails of the marker point profile. Finally, the resulting markerpoint profile is evaluated to determine its first moment, whichrepresents the marker point's location. Alternatively, the profile canbe processed by fitting a theoretical marker point profile to the markerpoint profile using a positional fit parameter. This method determinesthe position of marker points in a projection image of a computertomography scan with resolution significantly higher than that of theprojection image pixel size. Furthermore, the method performs markeridentification or localization in the presence of significantinhomogeneous background and adjacent structures. One method includesobtaining an image including a marker point, deriving a backgroundcorrected marker point profile, deriving a background and baselinecorrected marker point profile, and defining a sub-pixel position of themarker point as a center of a peak of the background corrected markerpoint profile.

Embodiments also provide systems and methods for rapidly identifying andassigning an approximate 3D position to a marker associated with markerpoints on the images. One method includes taking the observed markerpoint's positional changes from projection image to projection image,and fitting the positional changes to theoretical projectiontrajectories defined by using different possible values of the φ, Z, andR cylindrical coordinates. This step includes defining a set ofpotential values for one or more of the cylindrical coordinatedimensions and selecting one of the values from the set. The method thendetermines the marker points in each projection image that fall within arange of distances from the trajectory defined by the selected valuefrom the set. A metric is then assigned to each marker point falling inthe range by finding the mean squared distance of each marker point fromthe boundaries of the trajectory range. The metric is summed for allmarker points within the range to derive a single metric associated withthe particular trial coordinate value set. These steps are repeated foreach value in the set of potential coordinate values and the resultingvalues of the metric (as a function of coordinate values) are evaluatedto determine function peaks, which are each assumed to represent acluster of one or more markers. These coordinate values for the one ormore cylindrical coordinates associated with the peaks are assigned tothe corresponding derived physical markers. These steps are repeated forthe remaining cylindrical coordinates until the markers have beenassigned all three coordinates. Because of the nature of the metricfitting, it is possible to find the average marker location in thepresence of noise and patient movement and extract and separate multiplemarkers that are positioned close to each other. One method includesobtaining a sequence of images, each image representing an angle ofrotation of the scanner and having a first value and a second value. Themethod also includes analyzing the behavior of the first value of themarker point positions through the sequence of images to determine theradial distance of the marker and the angular position of the marker.The method also includes analyzing the behavior of the second value ofthe marker point positions through the sequence of images to determinethe axial position of the marker.

Embodiments also provide systems and methods for assigning each markerpoint in each projection image to the single physical marker most likelyto be associated with that marker point through a process called markertracking One method includes obtaining a sequence of images with aplurality of marker points and calculating an expected U and V positionand range for the physical marker in each image. The method alsoincludes selecting a physical marker and creating a list of candidatepoints made up of marker points that fall within a range of expectedmarker point locations in each projection image theoretically determinedfor that physical marker. The method further includes processing thelist by subtracting the theoretical trajectory of the physical markerfrom the actual image location of each marker point in the list toassign an error from theoretical to each marker point. The method thentraverses the list, projection-by-projection, and narrows the list tozero or one “good points” per projection image based on the first timederivatives between the errors of consecutive points in the candidatelist in both the U and V directions and other criteria.

Other aspects of the invention will become apparent by consideration ofthe detailed description and accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 schematically illustrates a cone-beam computed tomography system.

FIG. 2 schematically illustrates the computer of FIG. 1.

FIG. 3 illustrates an image including numerous marker points.

FIG. 4 is a flow chart illustrating a marker point identificationmethod.

FIG. 5 is a flow chart illustrating a background reduction process ofthe marker point identification method of FIG. 4.

FIG. 5A is a flow chart illustrating application of a 1-D cumulativemean filter.

FIG. 6 is a graphical illustration of a process used to rapidlymean-filter the image in one dimension that is part of the backgroundreduction process of FIG. 5.

FIG. 7 illustrates the image of FIG. 3 after applying the backgroundreduction process of FIG. 5.

FIG. 8 is a flow chart illustrating a marker point extraction process ofthe marker point identification method of FIG. 4.

FIGS. 9 and 10 illustrate the image of FIG. 7 after performing the stepsof the marker point extraction process of FIG. 8.

FIG. 11 is a flow chart illustrating a marker point sub-pixellocalization method.

FIG. 12 is a flow chart illustrating a marker point profile derivationprocess of the marker point sub-pixel localization method of FIG.

FIG. 13 is a graphical illustration of example pixel columns surroundinga marker point that are analyzed in the marker point profile derivationprocess of FIG. 12.

FIG. 14 is a graph illustrating example values of a marker point profileand a baseline profile used in the marker point profile derivationprocess of FIG. 12.

FIG. 15 is a flow chart illustrating a marker point positiondetermination process of the marker point sub-pixel localization methodof FIG. 11.

FIG. 16 is a graph illustrating a sub-pixel marker point center locationof an example marker point profile curve using the integral markercenter determination process of FIG. 15.

FIG. 17 is a graph illustrating a sub-pixel marker point center locationof an example best-fit Gaussian curve to a marker point profile curveusing the curve fit marker point center determination process of FIG.15.

FIG. 18 is a flow chart illustrating a physical marker three-dimensionallocalization method.

FIGS. 19A-B is a flow chart illustrating a process of analyzing a markerpoint's behavior in a U dimension.

FIGS. 19C-D is a flow chart illustrating a process of analyzing a markerpoint's behavior in a V dimension.

FIG. 20 is a graph illustrating an example of the process of fitting atrial value of φ to multiple U positions of a number of marker pointsthrough a sequence of projection images in the presence of motion in theU direction.

FIG. 21 is a graph illustrating an example one-dimensional plot of thecumulative fitting metrics for trial R and φ value pairs as a function φat a fixed value of R.

FIG. 22 is a graph illustrating an example of the process of fittingtrial values of Z to the V trajectories of three markers through asequence of projection images.

FIG. 23 is a graph illustrating an example one-dimensional plot of thecumulative fitting metric as a function Z.

FIGS. 24A-E is a flow chart illustrating a marker point to physicalmarker mapping method.

FIG. 25 is a graph illustrating an expected U range for an examplephysical marker and example U positions of all the marker pointsidentified through the sequence of projection images.

FIG. 26 is a graph illustrating the marker points illustrated in FIG. 25that have U positions within the U range illustrated in FIG. 25.

FIG. 27 is a graphical illustration of a portion of an example candidatemarker point list.

FIG. 28 is a graphical illustration of a flattening process of anexample candidate marker point list.

FIG. 29 is a graphical illustration of a synchronization process of anexample candidate marker point list.

FIG. 30 is a graphical illustration of processing a “next” projectionimage when the “next” projection image includes two or more candidatemarker points.

FIG. 31 is a graphical illustration of processing a “next” projectionimage when the “next” projection image includes a single candidatemarker point.

FIG. 32 is a graphical illustration of processing a “next” projectionimage when the “next” projection image includes no candidate markerpoints.

FIG. 33 is a graphical illustration of processing a “next” projectionimage when the good point list includes a good point and a suspectpoint.

DETAILED DESCRIPTION

Before any embodiments of the invention are explained in detail, it isto be understood that the invention is not limited in its application tothe details of construction and the arrangement of components set forthin the following description or illustrated in the following drawings.The invention is capable of other embodiments and of being practiced orof being carried out in various ways.

It should also be noted that a plurality of hardware and software baseddevices, as well as a plurality of different structural components, maybe utilized to implement the invention. Furthermore, and as described insubsequent paragraphs, the specific configurations illustrated in thedrawings are intended to exemplify embodiments of the invention.Alternative configurations are possible.

X-rays are a form of radiation similar to light rays but are rays thatcan penetrate the human body. When x-rays pass through a human body,however, the intensity of the x-ray decreases. The amount of decrease inintensity is related to the density of the tissues through which thex-rays pass. This characteristic of x-rays is used to produce images ofstructures inside the human body and other subjects.

An x-ray system traditionally includes an x-ray tube that generates adiverging beam of x-rays. A person or subject is placed between thex-ray tube and a sheet of film that is sensitive to the intensity of thex-ray striking it at each location. The x-ray tube is turned on and off,the film is developed, and the resulting image presents a “shadow image”of the structures in the path of the x-ray with denser structuresappearing whiter in the image.

Modern x-ray processes use one or more detectors rather than film. Thedetectors electronically detect and quantify the amount of (i.e., theintensity of) x-rays reaching the detector. Using the detectors, x-raycomputed tomography (“CT”) systems were created that rotate an x-raysource around a patient and electronically detect the resulting x-rayswith a single-wide strip of detector elements on the side of the patientopposite from the x-ray source. The data from the detectors is collectedfor all the different x-ray positions and then combined in a processknown as reconstruction. The combined images represent a single,finite-width “slice” through the patient where the image's intensity ateach point represents the x-ray density of the tissue at a particularphysical location. The combination of the x-ray source, the detectors,and the mechanical structure that allows rotation of the x-ray source isknown as the “CT gantry.”

Multiple slices can be acquired by repeating the process while movingthe patient, the x-ray source, or both between image acquisitions orscans. For example, moving the table supporting the patient while alsorotating the x-ray source produces a “helix” instead of a slice of data.In addition, increasing the size or width of the detector strip or ringfrom a signal row of detectors to multiple rows (e.g., up to 256 rows)allows for more rapid acquisition of more data. Furthermore, replacingthe detector strip with a larger two-dimensional detector acquires anentire detector panel image at each x-ray source position rather thanjust a single strip of data. The collections of these images, which cannumber 600 or more, are known as projection images. Each projectionimage represents an x-ray snap-shot of the patient from a differentperspective or angle as the x-ray source and the detector are rotated insynchronization around the patient. Because of the cone-shaped x-raybeam needed to cover the two-dimensional detector, this type of CTimaging is known as cone-beam (“CB”) CT imaging. FIG. 1 schematicallyillustrates a CB CT imaging system 10 according to one embodiment of theinvention.

The CB CT imaging system 10 includes a scanner 12 and a computer 14. Thescanner 12 includes a CT gantry 13 which includes an x-ray source 16, adetector 18, and a rotating carrier 20. The x-ray source 16 and thedetector 18 are aligned across from each other on the rotating carrier20, which moves the x-ray source 16 and the detector 18 around a patient22. The patient 22 is supported by a seat 24. The imaging system 10illustrated in FIG. 1 is a dental imaging system. Therefore, the patient22 sits in the seat 24 and places his or her chin into a rest 26. Therest 26 holds the patient's head relatively still while the gantry 13 isrotated to complete a scan of the patient's head.

The scanner 12 outputs image data to the computer 14. The image datarepresents the intensity levels of x-rays detected by the detector 18during the scan. The computer 14 is connected to a console 30, whichincludes a display 32 and one or more input and/or output devices, suchas a keyboard 34. A user uses the console 30 to interact with thecomputer 14. For example, a user can use the console 30 to requestimages or other data from the computer 14. The computer 14 provides therequested information to the console 30, and the console 30 displays theinformation on the display 32, prints the information to a printer (notshown), and/or saves the information to a computer-readable memorymodule (not shown).

FIG. 2 schematically illustrates the computer 14 of FIG. 1 according toone embodiment of the invention. The computer 14 includes aninput/output interface 40, an electronic processing unit (“EPU”) 42, andone or more memory modules, such as a random access memory (“RAM”)module 44 and a read-only memory (“ROM”) module 46. The input/outputinterface 40 receives image data from the scanner 12 and provides theimage data to the EPU 42. In some embodiments, the input/outputinterface 40 stores the image data to a memory module, such as the RAMmodule 44, and the EPU 42 obtains the image data from the memory modulefor processing. The input/output interface 40 may also transmit data tothe scanner 12, such as calibration and operational data orinstructions.

The EPU 42 receives the image data and processes the information byexecuting one or more applications or modules. In some embodiments, theapplications or modules are stored in memory, such as the ROM module 46.In other embodiments, the applications or modules may be stored on ahard disk storage unit and transferred to RAM for execution. As shown inFIG. 2, the ROM module 46 stores a marker point identification module50, a marker point sub-pixel localization module 52, a physical markerthree-dimensional (“3D”) localization module 54, a marker point mappingmodule 56, a patient motion identification module 58, and an imagereconstruction module 59. The EPU 42 retrieves and executes the markerpoint identification module 50 to identify marker points in the receivedimage data (see FIGS. 4-10). The EPU 42 retrieves and executes themarker point sub-pixel localization module 52 to provide accuratesub-pixel localization of marker point locations (see FIGS. 11-17). TheEPU 42 retrieves and executes the physical marker 3D localization module54 to determine three-dimensional locations of the physical markersidentified through their marker points in the projection images (seeFIGS. 18-23). The EPU 42 retrieves and executes the marker point mappingmodule 56 to map marker points identified in the projection images toparticular physical markers (see FIGS. 24-33). The EPU 42 retrieves andexecutes the patient motion identification module 58 to derive patientmotion information. The EPU 42 retrieves and executes the imagereconstruction module 59 to combine the projection images with themotion information and produce a set of motion-corrected axial images.

Complete gantry rotation takes approximately eight to 40 seconds. Duringthis time, a patient may move, which causes blurring of the resultingimages. Typical image resolution is on the order of 0.25 millimeter.Therefore, patient motion of this same order often causes image blurringand extensive patient movement can make the resulting imagesunacceptable for their intended clinical purpose.

The computer 14 (see FIG. 2) can correct the images for patient movementby tracking the movement of rigid objects in the resulting images. Forexample, in ideal conditions, where there is no patient movement, imagedrigid objects change location in the two-dimensions of the projectionimages in a well-defined way as the gantry rotates around the patient.Deviation between the expected locations of objects in the image iscaused by patient movement. By measuring a well-defined object'sdeviations from its expected locations, the amount of patient motion canbe measured and corrected. In particular, if at least three objects arepresent in an image, the measured deviations of the objects from theirexpected locations can be combined to determine a six-parameter patientposition error value (e.g., a patient motion vector), which can beapplied to the images to correct for the patient movement.

To ensure that the desired number of well-defined, rigid objects arepresent in the images, one or more fiducial markers are placed on apatient. The markers typically consist of lead or steel BBs, which aredense and prevent or limit x-rays from passing through. The markers maybe made from other materials and constructed in other forms or shapesthat are visible in a high proportion of projection images generatedduring a scan. Each marker or multiple markers may be positioned betweenlayers of adhesive and the adhesive can be applied to a patient toensure that the markers do not move during the procedure.

The markers are placed on a patient such that each field of view orimage created by the CB CT imaging system includes at least three markerpoints. For example, seven to nine markers may be placed on the patientto ensure that at least three marker points are within each imagecreated by the CB CT imaging system, to decrease the positionmeasurement noise, and to allow statistical combination of the results.In some embodiments, the markers are evenly spaced to have a maximalspatial distribution about the patient and to avoid interfering withnormal image interpretation.

Identifying Markers in an Image

After the markers are placed on the patient, a scan is made of thepatient, such as the patient's head, and the resulting image datarepresenting a sequence of projection images is transmitted to thecomputer 14 (see FIG. 2). FIG. 3 illustrates an example projection image100 that includes eight marker points 102. As an initial step incorrecting the generated projection images for patient movement, thecomputer 14 processes the generated projection images to identify themarker points in each image. FIG. 4 illustrates a marker pointidentification method 60 according to one embodiment of the invention.The marker point identification method 60 is performed by the EPU 42 ofthe computer 14 when the EPU 42 executes the marker identificationmodule 50.

As shown in FIG. 4, the first step of the method 60 generally includesreducing the intensity of the background of each projection image (step62). This background suppression step minimizes the background portionof the image. The background portion of the image includes everything inthe image other than the marker points. Because the fiducial markershave a well-defined shaped (i.e., a spherical or circular shape), step62 minimizes everything in the image that does not have thiswell-defined shape.

In some embodiments, background suppression is accomplished by applyinga linear high-pass filter to each image. In order to obtain acceptablespeed, three optimizations are performed. First, the high-pass operationis accomplished by combining the original image with a low-pass filteredimage. Secondly, instead of performing the filter simultaneously in thetwo image directions, a one dimensional high-pass operation issequentially applied to the image in two or more dimensions (i.e., at0°, 45°, 90°, and 135°). Finally, the low-pass filter operationaccomplished with a simple mean filter implemented with a cumulatedpixel approach described below.

A one-dimensional high-pass filter operation will tend to accentuateareas of rapid change in the direction of the filter. A sequence of suchoperations performed in different directions will tend to accentuatecircular objects over all other shapes. Since the marker points producedby the markers are circular, this sequence of filtering will accentuatethe marker points and attenuate any surrounding structures.

FIG. 5 illustrates the process or steps involved in the reduction ofbackground step 62. As shown in FIG. 5, the multi-dimensional high-passfiltering operation is executed by first performing simpleone-dimensional, cumulative mean, low-pass image filtering on theoriginal image A (step 69) to create filtered image B. Then, thefiltered image B is combined with the original image A to create image C(step 70) as described below. These two steps are repeated for eachfiltering direction using the output image from the previous filteringstep as the input to the next filtering step (steps 71-73). The processshown in FIG. 5 allows for rapid, multi-dimensional, high-pass filteringto reduce the background of the image. FIG. 6 illustrates an example ofthese steps involved in one-dimensional low-pass filtering in onedirection.

Rather than using a more typical kernel-based filtering approach, thefiltering of steps 69 and 71 is accomplished much more rapidly by usinga cumulative mean filter (depicted in FIG. 5A). This rapid filtering isaccomplished by creating a new image (image B1) where each pixel of thenew image is set to the cumulative sum of all pixel values from theoriginal image A from the particular pixel down to the first pixel inthe direction of the filter (step 74). For example, as shown in FIG. 6,when the computer 14 initially obtains an original image (image A) fromthe scanner 12, the method considers one row of image A that includesseven pixels (where the processing direction is horizontal), then acumulative image (image B1) is created by setting each pixel in this rowto the sum of the pixel values to the pixel's left plus this pixel'svalue.

Once the cumulative image is generated at step 74, image B1 is shiftedby an amount equal to the low-pass filtering factor (e.g., 2) to createa shifted image (image C1) (step 75). Image edges are handled by simplyzeroing shifted edge pixels. For example, as shown in FIG. 6, each pixelin the cumulative image (image B1) is shifted two positions to the leftto create the shifted image (image C1). Because there are no pixels twopositions left of the first two pixels (i.e., the edge pixels), thesepixels are set to 0.

The shifted image (image C1) is then subtracted from the cumulativeimage (image B1) to create a difference image (image D1) (step 76). Asshown in FIG. 6, to subtract the images, each pixel value in the shiftedimage (image C1) is subtracted from the corresponding pixel value in theun-shifted cumulative image (image B1). The final step is to normalizethe difference image (image D1) by dividing each pixel in it by thefilter factor (shift value) to produce the low-pass filtered image B(step 77).

The computation time for steps 74-77 is independent of the kernel sizeof the low-pass filter. The computations in these steps involve roughly2×N additions and subtractions instead of (kernel size)×N operations,which are typically needed for a normal kernel-based filter. For animage of 512×512 pixels and a typical smoothing factor of 15, thisrepresents a decrease in computation time versus a traditional 1-Dkernel approach of 7.5:1. Compared to a 2-D kernel, (versus two 1-Dcumulative mean filter operations) for a filter factor of 15, the speedimprovement is 56:1

After the one-dimensional low-pass filtered image is created, ahigh-pass filtered image (image C) is created by combining the low-passimage (image B) with the image used to generate it (image A) (step 73).This could be done by a simple difference operation between the twoimages. However, a division operation is used for the processing in thefirst direction filtering step 70. Division is used to compensate forthe exponential nature of x-ray attenuation and normalized image (imageB) to create image C (step 70). Dividing the difference image iseffective because x-ray attenuation occurs as a negative exponent of thematerial density. Therefore, because division is the equivalent ofsubtracting logs followed by exponentiation, dividing the differenceimage avoids the computationally expensive log and exponentiationoperations.

The process of applying the one-dimensional, cumulative-mean, low-passfilter and combining the input and output images is repeated in each ofthe successive dimensions (step 71-73). However, after the firstdimension is processed, a subtraction may be used instead of division atstep 72 because correction for the exponential nature of x-rayattenuation is completed by the division performed when the firstdimension is processed (step 70). Typically, four passes at 0, 45, 90,and 135 degrees are sufficient to obtain proper isolation ofspherical-shaped marker points from a wide range of backgrounds. FIG. 7illustrates the image 100 of FIG. 3 after applying the backgroundreduction method of FIG. 5. As shown in FIG. 7, the majority of thebackground (i.e., the portions of the image other than the marker points102) has been reduced or darkened, which allows the marker points 102 tobe more easily identified.

Returning to FIG. 4, once the background is reduced through the abovesteps, the marker points are extracted from the high-pass filtered image(step 80). FIG. 8 illustrates the process or steps of extracting themarker points from the high-pass image. As shown in FIG. 8, the processstarts by defining the local maxima in each high-pass image (step 82).This step includes identifying each pixel in each image in which all 8of the pixel's neighboring pixels have values smaller than it. Next, theidentified local maxima are used as seed points to grow regions (step84). This process combines the local maxima contained within aparticular region into a single candidate object or region. In someembodiments, the size or extension of the generated regions is limitedby an empirically-specified lower value threshold. FIG. 9 illustratesthe image 100 of FIG. 7 after region growing the local maxima (steps 82and 84), which creates candidate regions 104.

At step 86, a shape criterion is applied to each candidate region toidentify which regions are most likely to represent marker points. Insome embodiments, this process is performed by taking the ratio between(1) the pixel area of the candidate region (a) and (2) an area definedby the square of the mean radius of the candidate region times pi ((meanradius of a)²*π). The value of this ratio is used to determine whetherthe candidate region is the proper shape to be a marker point (e.g.,whether the candidate region has a circular shape).

Finally, at step 88, the candidate regions are further processed toreject unlikely candidates by comparing the location of candidateregions in adjacent projection images and applying a nearness criterion.For example, if a candidate region moves too much in its positionbetween sequential projection images, it is unlikely that the candidateregion represents a marker point. Similarly, if two candidate regionsare too close together, it is unlikely that the regions represent markerpoints. FIG. 10 illustrates the image 100 of FIG. 9 after applying theshape criterion to the candidate regions 104 and rejecting unlikelycandidate regions 104 (steps 86 and 88). As shown in FIG. 10, afterthese steps only candidate regions 104 corresponding to the eight markerpoints 102 (see FIGS. 7 and 8) remain.

In some embodiments, applying a two-dimensional high-pass filter to animage includes first applying a two-dimensional low-pass filter to theimage to create a low-pass image. Thereafter, the original image iscombined with the low-pass image.

In some embodiments, the marker identification method 60 can beperformed in times on the order of less than 50 milliseconds per image.In addition, when spherical markers are used, the marker pointidentification method 60 exhibits a high selective rate for markers anda high rejection rate for non-marker associated structures. For example,in some embodiments, the marker point identification method 60 can havea 100% sensitivity and a better than 40:1 specificity (less than 2.5%false detections) in real-world situations.

Determining a Sub-Pixel Center Point of a Marker

Even after the marker points are identified in the images using themarker point identification method 60, the marker points may need to befurther or more accurately defined before they can be used to performmotion correction. In particular, it may be necessary to know theeffective 2-dimensional position of the marker points to a resolutionmuch better than the pixel resolution (sub-pixel) provided by the markerpoint identification method 60. A resolution as small a ⅛ of the pixelresolution may be needed. Again, this must be done in the presence ofgreatly varying background both within the images and from image toimage. FIG. 11 illustrates a marker point sub-pixel localization method110 according to one embodiment of the invention. The marker pointsub-pixel localization method 110 is performed by the EPU 42 of thecomputer 14 when the EPU 42 executes the marker point sub-pixellocalization module 52. In some embodiments, the marker point sub-pixellocalization method 110 uses marker points identified by the markerpoint identification method 60 as a starting point to determine accuratesub-pixel marker locations as described below.

As shown in FIG. 11, the marker point sub-pixel localization method 110generally includes a marker point profile derivation process (step 112)and a marker point center determination process (step 114), which arerepeated for each of the two dimensions of the projection image (step116). As described in more detail below (see FIGS. 12-14), the markerpoint profile derivation process 112 applies a sequence of steps toremove the majority of the background and generate multi-pixel markerpoint profiles in each desired dimension. In the marker point centerdetermination process 114, the generated marker point profiles arecompared to the expected shape of the marker and one-dimensionalsub-pixel marker point center locations are derived (see FIGS. 15-16).

FIG. 12 illustrates the marker point profile derivation process 112according to one embodiment of the invention. Consider the operation ofderiving the vertical position of the marker point where a sequence ofpixels in the vertical direction is referred to as a column and asequence of pixels in the horizontal direction is referred to as a row.The process for deriving the horizontal position of the marker point isanalogous but with all references to rows replaced by columns andreferences to columns by rows. A row of pixels and column of pixels mayeach be referred to as a pixel line.

The first step of the process involves deriving a vertical backgroundprofile by averaging the pixel values over a range of background columnssurrounding the marker (step 120). FIG. 13 illustrates example columnsof pixels 121 surrounding and centered on a marker point 122. Step 120averages the pixel values of background columns 121 a surrounding themarker point 122 and excludes the pixel values of the marker columns 121b centered over the marker point 122. Averaging the background columns121 a includes summing the pixel values across each row of the columns121 a and dividing the sum by the number of columns. Averaging thebackground columns 121 a generates a single column of background values,which is a background profile for the marker point 122.

Next, a marker point profile is derived by averaging the marker columns121 b, which are the columns centered over the marker point 122 (step124). Again, averaging the marker columns 121 b includes summing thepixel values across each row of the columns 121 b and dividing the sumby the number of columns. Averaging the marker columns 121 b generates asingle column of marker values, which is a marker point profile for themarker point 122.

At step 126, the background profile is divided by the marker pointprofile on a row-by-row basis. For example, the value in each row of thecolumn of averages making up the background profile is divided by thevalue in the corresponding row of the column of averages making up themarker point profile. Dividing the row values in this order ensures thatimage inversion occurs (e.g., dark areas become light and light areasbecome dark). In addition, step 126 uses division (e.g., instead ofsubtraction) to compensate for the negative exponential nature of x-rayattenuation (since, as noted, division is mathematically equivalent totaking the exponent of the difference of the logs of the two profiles).

As an alternative embodiment, the log of the pixel values can becomputed before finding the marker point profile and background profile.Then background compensation of the marker point profile is accomplishedby subtracting the two profiles on an element-by-element basis andtaking the exponent of each element. This approach is slightly moreprecise than the above because it more accurately deals with imagesvariations over the averaged rows, but it takes more time than thesimple average and divide approach.

At this point, a reasonable baseline-corrected profile of themarker-point has been obtained. However, there may still be someresidual offset of tails of the profile from zero, particularly if thetails are asymmetric. The profile is brought to zero by the followingprocedure. An average baseline offset value is determined by averagingpixel values over columns in the tail distant from a peak of the markerpoint (e.g., a pixel within the marker point having the highest value),such as over columns where the pixel values are no longer decreasing asthey move away from the marker point center) (step 128). Care is takento exclude from this average pixels in the tails associated with anystructures adjacent to this marker point such as other markers or otherstructures in the image. FIG. 14 graphically illustrates the values ofan example marker point profile 129 and an example baseline profile 130and the results after fully baseline-correcting the marker point profile133.

Returning to FIG. 11, after the baseline-corrected marker point profilesare generated, the second step of the localization method 110 includesthe marker point position determination process 114. FIG. 15 illustratesthis process in more detail. As shown in FIG. 15, there are at least twoalternatives or options for performing this process. A first optionincludes determining the integral of the marker point profile (step 140)and defining the marker point position as the location along theintegral representing half of the integral (step 142). This marker pointposition is effectively the first moment of the integral, or theposition where that part of the integral to the left of this pointequals that part of the integral to the right of this point. FIG. 16illustrates an example sub-pixel marker point position location 144 ofan example marker point profile curve 146.

Returning to FIG. 15, the second option for determining a marker pointposition fits an appropriate curve function to marker point profileusing least-squares or other fitting (step 146) and defines the markerpoint position as the parameter representing the center of the curvefunction (step 148). FIG. 17 illustrates an example sub-pixel markerpoint position location 149 of a best-fit Gaussian curve 150 fitted toan example marker point profile curve 151.

The curve function used in step 146 may be a Gaussian curve or anexpected curve of a marker point, which takes into consideration issuesrelated to imaging physics, such as partial volume effects andexponential attenuation characteristics. In general, the curve functionused in step 146 may be any curve where one of the curve definitionparameters represents a curve center. The curve function may also handlenon-constant baseline profiles. In some embodiments, if a curve functionis selected that handles non-constant baseline profiles, steps 128 and132 in the marker point profile derivation step 112 can be eliminated.Furthermore, the curve function may be dynamically modified for eachmarker, image, or set of images based on measurable imagecharacteristics, such as known marker dimensions, pixel size, andapproximate sub-pixel location of a marker in other dimensions.

Determining 3D Cylindrical Coordinates for a Marker

Even after marker points and the marker point centers are identified inthe images (e.g., using the marker point identification method 60 andthe marker point sub-pixel localization method 110), the accumulation ofmarker points needs to be classified such that each marker point isproperly associated with the physical marker responsible for generatingthe point. In order to do this, the first step is to uniquely identifyeach physical marker point. The physical marker points are each uniquelyidentified through their three-dimensional (3D) coordinates. Therefore,the first step in associating marker points with physical markers is toidentify the physical markers and assign each its 3D coordinates. Thisprocess, however, can be complicated by patient movement and errors inthe determination of marker point position on the projection images. Inaddition, because there are multiple markers, some marker points may besuperimposed over other marker points in some projection images.

FIG. 18 illustrates a physical marker 3D localization method 160 foridentifying each of the physical markers and assigning each physicalmarker a position or location in 3D space. The physical marker 3Dlocalization method 160 is performed by the EPU 42 of the computer 14when the EPU 42 executes the physical marker 3D localization module 54.The physical marker 3D localization method 160 generally tracks theposition of the points produced by the marker through a sequence ofimages generated as the CT gantry rotates to map the marker to aposition in 3D space.

Existing tracking technologies use multiple sensors and multiple fixedrelationships between multiple detected sites to track the position of apatient in 3D space. The tracking technology or method applied in themethod 160, however, uses isolated single point detection sites.Therefore, rather than comparing the characteristics of multiple sitescollected at the same time, the tracking method collects characteristicsof single sites over time and matches the observed positional changeover time to a 3D position that would produce the same positional changeover time. This tracking method does not require any extra equipmentother than the CT imaging gantry and is tolerant of a relatively largeamount of image noise and patient movement.

As shown in FIG. 18, the method 160 starts by obtaining the position ofeach marker point (e.g., the center of each marker point determined bythe marker sub-pixel localization method 110) on each projection imagein two image dimensions, labeled “U” and “V” (step 162). U is adimension tangential to the circle of rotation of the CT gantry and V isa dimension parallel to the axis of rotation of the CT gantry. Thesedimensional positions of a marker point may be determined by the markeridentification point method 60, the marker point sub-pixel localizationmethod 110, or combinations thereof

If a marker has a fixed position in 3D space, its U and V positions (orvalues) on each projection image are described by the followingformulas:

$U = \frac{{DSD}*R*{\sin\left( {\theta - \phi} \right)}}{{DSO} - {R*{\cos\left( {\theta - \phi} \right)}}}$$V = \frac{{DSD}*Z}{{DSO} - {R*{\cos\left( {\theta - \phi} \right)}}}$where the 3D cylindrical coordinates (R=radial distance, φ=angularposition, and Z=height or axial position) represent the marker's 3Dposition, θ is the rotational angle of the CT gantry from its startingposition, DSD is the distance from the x-ray source to the detector, andDSO is the distance from the x-ray source to the object or patient.

To determine a marker's 3D position, the marker's cylindricalcoordinates or parameters (R, φ, and Z) must be determined using theabove U and V formulas and the marker's known U and V values. However,the U and V positions in a particular projection image may be differentfrom their correct values because of measurement error. Also, R, φ, andZ are functions of time (e.g., R(t), φ(t), and Z(t)) due to patientmotion.

Due to these difficulties, the method 160 finds “average” values for R,φ, and Z (step 164). The average values are those values that producethe minimum sum squared error of the marker's actual position in eachprojection image from the marker's expected position in each projectionimage as specified by the above U and V formulas. The following errorformula generally provides these average values:

${Err} = {\sum\limits_{\theta = 1}^{n}\;\left( {\left\lbrack {U_{meas} - \frac{{DSD}*R_{avg}*{\sin\left( {\theta - \phi_{avg}} \right)}}{{DSO} - {R_{avg}*{\cos\left( {\theta - \phi_{avg}} \right)}}}} \right\rbrack^{2} + \left\lbrack {V_{meas} - \frac{{DSD}*Z_{avg}}{{DSO} - {R_{avg}*{\cos\left( {\theta - \phi_{avg}} \right)}}}} \right\rbrack^{2}} \right)}$

However, the process of finding the desired values for R_(avg), Z_(avg),and φ_(avg) is simplified because the U dimension is independent of Z.As such, the process can be divided into two processes or steps: (1)analyzing the marker's behavior in the U dimension (see method 166 inFIGS. 19A-B) and (2) analyzing the marker's behavior in the V dimension(see method 168 in FIGS. 19C-D).

FIGS. 19A-B illustrates a method 166 of analyzing the marker's behaviorin the U dimension according to one embodiment of the invention. Asshown in FIG. 19A, the method 166 starts by selecting a number set of Rand φ value pairs representing a sub-sample of all possible value pairsfor the R and φ parameters (step 170). The size of the sub-sample setrepresents a compromise between execution speed and obtainable accuracyfor R and φ. In general, the R and φ values are chosen to uniformlysample the range of physically possible R and φ values. That is, R isbetween 0 and the maximum scanner radial field of view, and φ is between0° and 360°. Alternatively, a value of R may be assumed and only asub-sample of possible values of φ is selected. Therefore, although theremaining steps of the method 166 are described with reference to pairsof R and φ values, it should be understood that the steps can beperformed with an assumed value of R and varied values of φ.

As shown in FIG. 19A, after the sub-sample set is selected, an R and φvalue pair is selected from the sub-sample set to process (step 172).Next, using the selected R and φ value pair and the U formula describedabove, a “theoretical U” value is determined for each projection image(step 174). Because each projection image is associated with aparticular angle of rotation of the CT gantry, the θ value associatedwith each projection image is used in the U formula to determine each“theoretical U.” A U range is then defined around each “theoretical U”value (step 176). The U range is centered about the “theoretical U”value and has a magnitude in each direction from the center equal to themaximum expected U-variation of a marker from its average or expected Uposition caused by patient motion.

Next, the U positions of the marker points identified in each projectionimage are compared with the U range for each projection image. In step178, an identified marker point is selected. If the selected markerpoint identified in a projection image has a U value that falls withinthe projection image's U range (as determined in step 180), the markerpoint is identified as a marker point candidate and the method proceedsto step 182. Marker points that fall outside this range are neglected.It should be understood that more than one marker point may have a Uvalue within a U range for a particular projection image. In fact, ifthere are multiple adjacent physical markers positioned on a patient, itis likely that more than one marker identified in the projection imageswill have a U value within a U range for a particular projection image(i.e., a particular value of θ).

FIG. 20 illustrates the U positions of a number of markers through asequence of 300 projection images. As shown in FIG. 20, three trials orvariations of φ values were performed, while a corresponding R value wasassumed and held constant. In particular, a first trial used a φ valueof 60 degrees, a second trial used a φ value of 100 degrees, and a thirdtrial used a φ value of 140 degrees. The U ranges for each of thesetrials are shown in the graph with dashed lines. For each projectionimage (see x-axis), each marker point that has a U position within the Uranges will be marked as a candidate marker.

Next, a fitting metric is determined for each candidate marker point ineach projection image (step 182). In some embodiments, the fittingmetric is the square of the distance between the marker candidate's Uvalue and the closer boundary of the U range. For example, if a markerhas a U value of 20 and the U range is 0 to 30, the fitting metric wouldbe set to 10 (which is the distance between 20 and 30—the closer U rangeboundary) squared or 100. If the marker's U value is the same as the“theoretical U” value, the fitting metric is maximized. The fittingmetrics values from all candidate points are then summed to create acumulative fitting metric, which is assigned to the currently selected Rand φ value pair (step 184). As shown in step 186 of FIG. 19A, steps178-184 are repeated for all points in all projection images.

As shown in FIG. 19A, steps 172-186 are repeated for each R and φ valuepair in the sub-sample set (step 188) to assign a cumulative fittingmetric to each value pair. After all of the R and φ value pairs havebeen assigned cumulative fitting metrics, the method 166 evaluates thecumulative fitting metrics for all R and φ value pairs as a function ofR and φ (step 190 in FIG. 19B). The method 166 uses peak detection todetermine the peaks of R and φ in this function (step 192). The peaks ofthe function are the peak R and φ value pairs which each represent oneor more physical markers. The method 166 assigns each peak R and φ valuepair to the corresponding marker associated with that peak (step 194).FIG. 21 illustrates an example of a one-dimensional plot (e.g., ahistogram) of the cumulative fitting metrics for the R and φ value pairsas a function of φ (the value of R is assumed and held constant). Asshown in FIG. 21, the peaks of the plotted function represent φ valuesfor the physical markers.

After analyzing the marker points' behavior in the U dimension, althoughcandidate markers have been located and assigned R and φ values, it ispossible that more than one marker exists with the same R and φ values.To address this situation, the marker points' behavior in the Vdimension is analyzed. FIGS. 19C-D illustrates a method 168 of analyzingthe marker points' behavior in the V dimension according to oneembodiment of the invention. As shown in FIG. 19C, the method 168 startsby utilizing the peak R and φ value pair assigned to that markercandidate during method 166 (see step 194 in FIG. 19B). It uses thesevalues to identify candidate points to evaluate further in terms of Z(step 200). A sub-sample of all possible values of Z is then selected.(step 202). The number of values in the sub-sample is chosen as atrade-off between execution speed and desired Z-value resolution. Therange of the sub-sample is chosen to cover the range of physicallypossible values of Z (points that can at some point be within thescanning field of view).

A Z value is then selected from the sub-sample set (step 204). Using theR and φ value pair previously determined for this marker or markers, theselected Z value, and the V formula described above, a “theoretical V”value is determined for each projection image (step 206). Because eachprojection image is associated with a particular angle of rotation ofthe CT gantry, the θ value associated with each projection image is usedin the V formula to determine each theoretical V. A V range is thendefined around each theoretical V value (step 208). The V range iscentered about the theoretical V value and has a magnitude in eachdirection from the center equal to the maximum expected V-variation of amarker from its average or expected V position caused by patient motion.

Next, the V position of each marker point that is part of the candidatelist has its V value compared against the V range for this Z trial. Instep 210, a point in the projection image is selected. If the point hasa V value that falls within the projection image's V range (step 212),the marker is identified as a further restricted marker point candidate.To clarify, a point must fall within the range of U values associatedwith its R and φ values and within the range of V values associated withits R, φ, and Z values in order to be considered a candidate. Asdescribed above, more than one marker point identified in a particularprojection image may have V and U values within the current V and Uranges associated with that projection image. As shown in step 216 ofFIG. 19C, steps 210-215 are repeated for all points in all projectionimages.

FIG. 22 illustrates the V positions of three markers through a sequenceof 300 projection images. In FIG. 22, three trials or variations of Zvalues are shown. In particular, a first trial used a Z value of 70millimeters, a second trial used a Z value of 30 millimeters, and athird trial used a Z value of −15 millimeters. The V ranges for each ofthese trials are shown in the graph with dashed lines. For eachprojection image (see x-axis), each marker that has a V position withinthe V ranges will be marked as a marker candidate.

Returning to FIG. 19C, a fitting metric is determined for each markercandidate in each projection image (step 214). The fitting metric is thesquare of the distance between the marker candidate's V value and thecloser boundary of the V range. For example, if a marker has a V valueof 20 and the V range is 0 to 30, the fitting metric would be set to 100(which is 10—the distance between 20 and 30, the closer V rangeboundary—squared). The metric fittings from all the projection imagesare then summed to create a cumulative fitting metric, which is assignedto the currently selected Z value (step 215).

As shown in FIG. 19C, steps 204-218 are repeated for each Z in thesub-sample set (step 217) to assign a cumulative fitting metric to eachvalue. After all of the Z values in the sub-sample set have beenassigned cumulative fitting metrics, the method 168 evaluates therelationship between Z and the resulting cumulative fitting metric (step218 on FIG. 19D). The method 168 then uses peak detection to determinethe peaks of the relationship (step 220).

Each peak in the cumulative fitting metric—Z relationship is associatedwith a unique physical marker. The Z value of the peak represents the Zvalue of the corresponding physical marker. The R and φ values assignedare those previously found for the marker or markers. Each of thesemarkers now has known values for all three cylindrical coordinates (R,φ, Z). FIG. 23 illustrates an example of a one-dimensional plot (e.g., ahistogram) of the cumulative fitting metric as a function of Z. In FIG.23, there are three peaks corresponding to three separate physicalmarkers, each at the same value of R and φ.

The method 168 is (steps 200-222) is repeated (step 224) for each of theother R and φ pairs determined by method 166 to isolate all uniquephysical markers and determine Z values for each. The method 168terminates when all marker candidates are separated into unique physicalmarkers and each physical marker has been assigned all three cylindricalcoordinates. The assigned cylindrical coordinates represent the marker'saverage R, φ, and Z values, which represent the marker's 3D position.

The methods 166 and 168 described above with respect to FIGS. 19A-Drepresent a two-step process which uses two separate methods to assign3D cylindrical coordinates to a marker (i.e., a first method assigns theR and φ values and a second method assigns the Z value). A firstalternative implementation of the marker 3D localization method 160combines these methods into a single operation where a sub-sampled setof values of R, φ, and Z is selected and each combination of values isused to determine a single fitting metric based on the combined U and Verror. The result of this single process is a three-dimensional arraythat has peak values coincident with each isolated physical marker. Thisimplementation, however, may take more time than the separate methodimplementation. In particular, if R, φ, and Z have, respectively, m, n,and q sub-samples, the first alternative implementation requires m*n*qiterations versus (m*n)+q for the separated implementation.

A second alternative implementation of the marker 3D localization method160 can also be used to speed up the separated implementation. In thesecond alternative implementation, R, φ, and Z are found as a two-stepprocess as originally described, but the number of R values (m)processed in the first method 166 would be decreased by sub-sampling Rat a coarser frequency. This is reasonable because the fitting is lesssensitive to R values. Then, a third method can be added that holds φand Z fixed and uses a finer-sampled, sub-sample set of R and thefitting metric process described above to optimize R.

A third alternative implementation of the marker 3D localization method160 can also make the sub-sampling coarser for one parameter, twoparameters, or all three parameters. This implementation can alsoinclude final fitting methods for each marker candidate that uses afiner sampling rate for each parameter sub-sample set covering anarrower range. This implementation may increase both the speed of themethod 160 and improve the final accuracy of the marker positioncoordinates.

One general challenge for the marker 3D localization method 160 ishandling a situation where the markers are positioned very close to eachother but still allow a reasonable range of motion for each marker. Ifthe markers are separated by less than the maximum expected motion, themethod 160 might label two markers as if they were one marker and assignthis one marker a position between the actual positions of the twomarkers. To prevent this situation, the markers should be placed on thepatient with adequate physical separation. But the method 160 needs todeal with situations where placement guidelines are not followed. Themethod 160 can be modified to identify situations where a largepercentage of projection images have more than one marker within thecandidate point acceptance range or identify situations where acumulative metric peak is larger than expected for a single marker. Ifeither of these situations is detected, the method 160 can includeadditional processing to better isolate the two markers. One approachwould be, when overlapping markers are suspected, to repeat thelocalization process with a finer sub-sampling and ranges narrowed tocover only the suspect marker.

A second challenge is when marker motion causes two peaks in thecumulative fitting metric. In this case, a single marker will appear astwo. To address this possibility, the method 160 can divide the θ rangeinto sub-ranges, perform fitting within each sub-range separately, andthen merge the results. In addition, the method 160 could adaptivelychange the candidate point acceptance range for the three coordinateparameters, R, φ, and Z. For example, the method 160 could start with asmall acceptance range to ensure adequate isolation of closely adjacentmarkers and assuming little motion. If an adequate number of points arenot included in these ranges, presumably because of motion, the rangecould be increased until sufficient candidate marker points are found.

In general, the marker 3D localization method 160 localizes markers tocoordinate resolutions better than one millimeter in the presence ofpatient motion as great as two centimeters over the course of the scan.The method 160 also operates effectively in the presence of 24 or moremarkers and is tolerant of measurement noise. Finally, the method 160can process 24 markers over 300 projection images in less thanapproximately 5 seconds.

Mapping a Point to a Physical Marker in a Patient

Even after image marker points are identified and localized and thecorresponding physical markers identified and assigned 3-D coordinates,each marker point must be assigned to a specific physical marker placedon the patient. In a typical CB CT scan with 600 images and ninemarkers, there will be over 5,000 identified points in the projectionimages that must each be assigned to one of the nine physical markers.

Looking at the marker points identified in a single projection image,there is not (at least not in general) enough information to assign aparticular marker point identified in an image to a particular physicalmarker positioned on the patient. In the ideal case, if the markers arephysically motionless and their physical locations known exactly, uniqueU and V tracks of marker point position change from one projection imageto a successive image can be defined for any marker point identified inthe projection images. Only one physical marker will be consistent withthe point position track that a particular marker makes from oneprojection image to the next. Thus, two projection images aretheoretically sufficient to assign these two marker points to theircorresponding physical marker.

However, additional challenges can occur when either a marker's physicallocation in the projection images is not exactly known or if there is asignificant amount of patient motion (or unexpected gantry motion)during the scan. In these cases, a level of trajectory uncertainty isadded that complicates the proper assignment of a marker pointidentified in projection images to a physical marker placed on apatient.

FIGS. 24A-E illustrates a physical marker mapping method 250 forassigning marker points identified in projection images to specificphysical markers placed on a patient according to one embodiment of theinvention. Generally, the method 250 tracks marker point locations insuccessive projection images in the presence of patient motion andevaluates each marker point on a projection image-by-projection imagebasis to identify the marker points most likely to be associated with aparticular physical marker. The method 250 uses a number of differentcriteria, including the distances between the U and V measured for amarker point and the theoretical U and V associated with a physicalmarker and the time derivatives of these distances.

As shown in FIG. 24A, an initial step of the method 250 includesidentifying all physical markers associated with the scan and obtainingan approximate 3D location for each (step 252). This can be done usingphysical measurements, measurements from the reconstructed imagedataset, or direct analysis of the projection images. In someembodiments, the method 250 uses the physical marker localization method160 described above to obtain an approximate 3D location for eachphysical marker.

Next, the method 250 chooses one physical marker and calculates anexpected U position (U_(image)) and an expected V position (V_(image))for each projection image (step 254). The method 250 uses the followingequations that describe how a physical point in space is expected tomove from image to image based on the particular geometry of the CTgantry.

$U_{image} = \frac{{DSD}*R*{\sin\left( {\theta - \varphi} \right)}}{{DSO} - {R*{\cos\left( {\theta - \varphi} \right)}}}$$V_{image} = \frac{{DSD}*Z}{{DSO} - {R*{\cos\left( {\theta - \varphi} \right)}}}$

For example, assuming that the CT gantry rotates smoothly about a centerof rotation, R is the distance a point is from this center in a plane ofrotation, φ is the angular location of the marker, and Z is the distancethe marker is from the central plane of rotation (i.e., the plane thatincludes the point defining the location of the X-ray source). DSD isthe distance from the x-ray source to the detector plane, and DSO is thedistance from the x-ray source to the CT gantry's center of rotation.Each projection image is acquired at a gantry rotation angle of θ. Usingthese equations, the method 250 defines the expected U and V projectionimage positions of a physical marker with a 3D position of R, φ, and Zat each projection image angle, θ.

The method 250 then uses the expected U and V marker image positions todefine an expected U range and an expected V range of trajectories forthe physical marker (step 256). The ranges are large enough to includeall points associated with the physical marker, even in the presence ofpatient motion. For example, if the patient is expected to move as muchas 20 millimeters over the course of the image acquisition, the expectedU and V ranges for a particular physical marker would be:

${\frac{{DSD}*R*{\sin\left( {\theta - \varphi} \right)}}{{DSO} - {R*{\cos\left( {\theta - \varphi} \right)}}} - 10} < U_{image} < {\frac{{DSD}*R*{\sin\left( {\theta - \varphi} \right)}}{{DSO} - {R*{\cos\left( {\theta - \varphi} \right)}}} + 10}$and${\frac{{DSD}*Z}{{DSO} - {R*{\cos\left( {\theta - \varphi} \right)}}} - 10} < V_{image} < {\frac{{DSD}*Z}{{DSO} - {R*{\cos\left( {\theta - \varphi} \right)}}} - 10}$

The method 250 then creates a candidate point list for the physicalmarker (step 258) and adds each marker point identified in eachprojection image that falls within both the expected U range and theexpected V range (step 260). FIG. 25 is a graph illustrating an exampleexpected U range for a particular physical marker (between the dottedlines) and example U positions of marker points identified through asequence of projection images. FIG. 26 illustrates those markersillustrated in FIG. 25 that have U positions within the expected U rangeillustrated in FIG. 25. Any of these markers that also have a V positionwithin the expected V range will be added to the physical marker'scandidate point list. For example, FIG. 27 graphically illustrates aportion of a candidate point list. As shown in FIG. 27, the candidatepoint list includes each marker from each projection image that has a Uposition and a V position that falls within the expected U and V ranges.The final list of marker points associated with the physical marker willbe selected from this candidate point list. As shown in FIG. 27, thecandidate point list can include more than one marker within aparticular projection image (two candidate points in this example occurboth in projection 1 and project 6. This can occur when two markerscross each other's path through the sequence of projection images. Therecan also be more than one marker point in a candidate point list for aparticular projection image if two physical markers are positioned closeto each other on the patient.

As shown in FIGS. 24A and 24B, after creating and populating thecandidate point list (steps 258 and 260), the method 250 processes thecandidate point list by subtracting the “expected” marker locations ofthe physical marker from the actual marker point's locations for eachpoint included in the candidate point list for both the U and Vdimensions (step 262). This process flattens the marker's locations tocreate an “error from expected” candidate point list. The “error fromexpected” candidate point list, therefore, represents the amount ofdiscrepancy or error (e.g., caused by patient movement or measurementerror) between the expected position of the physical marker reflectingin the projection image and the actual marker point position. FIG. 28graphically illustrates this flattening process applied to an examplecandidate point list.

Once the “error from expected” candidate point list is created for aphysical marker, the process of point tracking occurs. The objective ofthis process is to proceed from one projection to the next and selectthe one point in the next projection most “likely” to be the correctpoint, that is, the point most likely to have been truly created by theassociated physical marker. A number of heuristics are applied in orderto accomplish this. The overall process generally has two components:synchronization and tracking In the synchronization component, the pointlist is processed, one projection at a time, until a particular point isfound to have a high likelihood of being associated with the physicalmarker. The tracking component starts with this point and continues,projection by projection, selecting the one (or zero) point (from thelist of zero to multiple points at that projection) most consistent witha track that would be expected to be formed by this marker. If noconsistent points are found for a set number of images, synchronizationis lost and will have to be re-obtained. The tracking process canproceed in either direction of increasing projection number ordecreasing projection number.

In the tracking process, at each new projection there are three possiblecases: zero points exist, one point exists, and more than one pointexists. Motion and measurement noise complicates the process in that a“correct point” might not be where it is expected to be. This is handledby implementing the concept of “suspect points.” If a best candidatepoint falls outside some criteria, it is labeled as suspect. This labelis removed and the point either included or rejected by looking at itsposition in the context of point positions in future projections. Thefollowing provides specifics for this process.

In some embodiments, the tracking component of method 250, steps266-352, may be performed in either the direction of increasingprojection number or that of decreasing projection number. In the lattercase, it should be understood that references to “next projection” and“previous projection” and other terms relating to sequencing arerelative to the direction of processing.

Returning to FIG. 24B, after processing the candidate point list at step262, the method 250 creates an empty good point list for the physicalmarker (step 264). As described in more detail below, the good pointlist is populated with marker points identified in the projection imagesthat have a high likelihood of being associated with the physicalmarker. The first step is synchronization (steps 266-270). The method250 examines the candidate list, one projection at a time (e.g.,starting from the first projection image) until a defined number ofsequential projection images, such as 3, are found, where each image hasa single candidate point in the candidate point list (step 266). Thepoint associated with the first projection image in the detectedsequence of projection images is labeled as the “reference good point”or reference point (step. 268). In addition, the method 250 adds thispoint to the good point list (step 270). FIG. 29 graphically representsthe synchronization process on an example candidate point list. As shownin FIG. 29, the candidate point list includes a sequence 271 of threesequential projection images, where each image in the sequence 271 hasonly one point in the candidate point list. In the synchronizationprocess, the first projection image in the detected sequence (labeled271 a) contains the single candidate point 271 b which is labeled as thereference good point (step 268) and added to the good point list (step270).

After synchronization has been obtained, the method 250 then moves tothe “next” projection image in the sequence (step 272). In someembodiments, the method 250 will use the projection image found in thesynchronization process and proceed first backward then forward fromthis point. The backward processing will proceed either until aprojection having a point in the good point list is reached or, in thecase where this is the first synchronization, to the first projection.The forward process will proceed until either the last projection isreached or until synchronization is lost. Again, when processingsequential projections, the definition of “next” will depend on whetherprocessing is forward or backward.

After moving to the next projection image, the method 250 determines howmany candidate points are in that projection image (step 274). In someembodiments, there are two categories of results for step 274. The nextprojection image can include one or more candidate points or zerocandidate points. FIG. 30 illustrates a situation where a nextprojection image 275 has two candidate points. FIG. 31 illustrates asituation where a next projection image 275 has only one candidatepoint, and FIG. 32 illustrates a situation where a next projection image275 has no candidate points.

Returning to FIGS. 24B and 24C, if the next projection image has one ormore candidate points, the method 250 determines the derivative of the Uand V positional error with respect to θ (projection angle) between eachcandidate point in the next projection image and the most recently-added“good point” in the good point list (step 276). This derivativerepresents, in effect, the flatness of the line connecting two points asin FIG. 27. A zero derivative means that the amount of error betweenactual and theoretical has not changed between two images.

Next, the method 250 selects the candidate point associated with thederivative closest to zero (step 278) and adds this point to the goodpoint list (step 280). If there is only a single candidate point, thispoint is chosen. Because the candidate point list was flattened toremove all effects except error (e.g., caused by patient movement ormeasurement noise), a large variance in the position error derivativebetween marker points in two adjacent projections may indicate (1) thata structure or object in an image has been erroneously defined as amarker or (2) that the marker point is actually associated with adifferent physical marker.

Therefore, two markers in the candidate point list that are associatedwith adjacent projection images and have a derivative close to zero arelikely to be adjacent markers associated with the same physical marker.For example, as shown in FIG. 30, the next projection image 275 includestwo candidate points 292 a, 292 b. The derivative between the reference“good point” 293 and the candidate point 292 a is closer to zero thanthe derivative between the reference good point 293 and the candidatepoint 292 b. Therefore, the candidate point 292 a will be selected atstep 278. As another example, as shown in FIG. 31, the next projectionimage 275 includes one candidate point 292 a. The error derivativebetween the reference good point 293 and the candidate point 292 a iscalculated and the candidate point 292 a is added to the good point listas a good point at step 280.

After selecting the candidate point associated with the derivativeclosest to zero, the method 250 determines if the candidate point'sderivative passes a derivative test (step 282). The derivative test mayinclude determining if the derivative has a value that differs from zeroby more than a predetermined amount, such as 0.02 mm. If it does not,the candidate point passes the derivative test and the method 250 labelsthe point as the reference good point (step 284). If the point does notpass the derivate test, the method 250 labels the point as a “suspectpoint” (step 286). This newly added point is used to process thecandidate points in the subsequent projection image.

The predetermined amount for the derivative test can be set and changedto modify the sensitivity of the method 250. For example, thepredetermined amount can be set to a larger value in order to allowfewer suspect points to be added to the good point list, which mayincrease the amount of marker points ultimately mapped to a particularphysical marker. In some embodiments, the derivative test can vary basedon the number of candidate points in the next projection image, whetherthe good point list includes a suspect point, etc. In some embodiments,points labeled as suspect may not be considered for future processingbut rather always considered bad and not added to the good point list.In these cases, the associated projection would have zero points.

Returning again to step 274 of FIG. 24B, if the next projection imagedoes not include any candidate points, the method 250 skips the nextprojection image and continues to use the most recently-added referencegood point (and the most recently-added suspect point if applicable) toprocess candidate points in the subsequent projection image (step 292).For example, as shown in FIG. 32, the next projection image 275 does notinclude any candidate points. Therefore, the method 250 skips the nextprojection image 275 and continues to use the most-recently addedreference good point 293 (i.e., from the current projection image 295)to process the candidate point(s) (e.g., a candidate point 296 a) in thesubsequent projection image 297.

A suspect point can exist either because measurement noise causesmisplacement of a particular point or because patient or unplannedgantry motion has caused actual movement in the physical marker. Forproper motion correction, it is desirable to reject points related tothe former and include points related to the latter. The handling ofsuspect points addresses these motion correction goals. Basically, theimplementation of schemes for handling suspect points is based on theconcept that a point that is suspect due to motion likely will havebehavior reflected in multiple adjacent projections. Thus, once asuspect point has been flagged, final classification depends onexamining future and possible additional past projection behavior.

As referred to above, if the method 250 determines in steps 288 or 290that there is a previous suspect point in the good point list (alsoreferred to as an old suspect point), and the next projection image hasone or more candidate points, special processing is required (see steps300 and 302 in FIG. 24D). If so, the error derivative is formed betweenthe suspect point and the most recently added point in the good pointlist (step 300 and 302). After forming the derivative associated withthe suspect point, the method 250 determines if this derivative passes aderivative test (steps 304 and 306). The derivative test may includedetermining if the derivative has a value that differs from zero by morethan a predetermined amount, such as 0.02 mm. If the derivative differsby less than the predetermined value, the suspect point is said to passthe derivative test. If the derivative differs by the same or greaterthan the predetermined value, the suspect point is said to fail thederivative test. If this suspect point passes the derivative test, thesuspect label is removed from the point (steps 308 and 310). If thesuspect point fails the derivative test, this point is removed from thegood point list (steps 312 and 314) since it has failed two derivativetests, one with the previous point and one with the succeeding point.

If the previous suspect point passes the derivative test and point mostrecently added to the good point list had failed its derivative test,this previous suspect point now becomes the reference good point in theprocessing of succeeding projections (step 316). In this case, it isassumed that motion had occurred and the previous suspect point nowrepresents the best indication of the new marker position.

For example, as shown in FIG. 33, the first derivative between thesuspect point 321 and the candidate point 323 of projection image 322 isclosest to zero and passes the derivative test. In addition, thecorresponding first derivative between the good point 320 and thecandidate point 323 fails the derivative test. Therefore, the suspectpoint 321 has its suspect label removed at step 310. The previouslylabeled suspect point, 321, will be used as the reference point forprocessing the next projection image (step 316).

There are a number of alternative embodiments of point tracking, manydealing with alternative ways of handling suspect points. In one, themethod 250 is configured to only accept or re-label suspect points asgood points if they have an error derivative with a candidate point thatis closer to zero than the corresponding derivative with the referencegood point for a predetermined number of sequential projection imagesgreater than one. For example, the method 250 may determine whether, forthree projection images in a row, the derivative between the suspectpoint and a candidate point is closer to zero than the correspondingderivative between the reference good point and the candidate point. Ifso, the method 250 re-labels the suspect point as a good point. Thenumber of required projection images for accepting a suspect may be setand modified to change the sensitivity of the method 250 (i.e., thenumber of suspect points that ultimately get accepted as good points).In another alternative embodiment, any suspect points could beimmediately removed from any further consideration. In such cases, theremight be more reliance on resynchronization. In another embodiment, asuspect point flag would replace the labeling of individual points assuspect. In this case, the processing performed with the suspect pointflag is set might be different depending on whether the candidate pointlist has a single point or more than one point.

If more than a particular number of projections, such as five, have nocandidate points as determined in step 294, the processing is declaredto have lost synchronization (step 299). In these cases, thesynchronization process described above is repeated to regainsynchronization.

As shown in FIGS. 24C and 24D, after processing the candidate point(s)in the next projection image, the method 250 defines the next projectionimage as the new current projection image (step 350) and repeats steps272-350 (i.e., those steps applicable) for the new current projectionimage. This process is repeated until all the projection images in onedirection from the initial current projection image have been processedas the “next projection image” (e.g., the first projection image isreached and processed as the next projection image) (step 352).

In some embodiments, the method 250 projection image steps 266-352 maybe completely performed twice, once going from the lowest projectionnumber to the highest and once going from the highest projection numberto the lowest with each direction generating a separate good point list(step 354). Once the point lists are generated, they are combinedprojection-by-projection into a new list with any points not in bothoriginal lists eliminated (step 355). This combination insures that anypoints included in the final good point list have a higher likelihood ofbeing correct. FIGS. 30-32 illustrate processing projection images in abackward or reverse direction. FIG. 33 illustrates processing projectionimages in a forward direction.

After the projection images have been processed in either direction, themethod 250 repeats steps 254-354 for each physical marker (step 356).For example, if nine physical markers were placed on the patient beforethe scan, the method 250 repeats steps 254-354 for each of the ninephysical markers. After steps 254-354 have been performed for eachphysical marker, the method 250 has created a good point list for eachphysical maker and each good point list includes “good points”representing those marker points identified in the projection imagesthat have the highest likelihood of being associated with a particularphysical marker (step 360). Each good point list can include only zeroor one marker point from each projection image that is associated with aparticular physical marker. The “good points” in the good point listscan be used to extract patient motion information. For example, theexpected locations of a physical marker reflected as marker pointsthrough the sequence projection images can be compared to the actuallocations of the markers associated with the physical marker through thesequence projection images to determine the amount of patient movement.

One general challenge is that if physical markers are positioned closeto each other (e.g., closer together than the range allowed for patientmovement), it may be difficult to properly assign markers to aparticular physical marker. This can be avoided by ensuring that thephysical markers are separated by a proper distance when they are placedon the patient. Nonetheless, the method 250 is designed to handle, atleast in part, situations where placement guidelines are not followed.If this situation occurs, it is possible that the tracking process mightincorrectly assign marker points, a process referred to as “trackjumping.” Once one marker point “jumps” to a wrong track, there may notbe a good basis for jumping back to the correct track, unless a sequenceof points in the new (and wrong) track are missing. While the method 250can deal with missing points on a track of markers, excessive missingpoints can lead to track jumping. This possibility can be minimized bychoosing a small range of candidate points for a physical marker'scandidate point list (e.g., because of lower expectations of the amountof patient movement).

In addition, processing the projection images in both directions (e.g.,forward and backward) can identify and eliminate track jumping. Forexample, by comparing the results of processing the projection images inthe forward direction to the results of processing the projection imagesin the backward direction, the method 250 can identify extended lengthdifferences, which are an indication of tracking problems. Another wayto identify tracking problems is to compare the good point lists fordifferent physical markers. If two good point lists include more than afew (e.g., one or two) of the same markers, the method 250 can detectand correct an incorrect tracking problem. However, with reasonablemarker point detection and physical markers that are separated byapproximately one centimeter or more, track jumping is usuallynon-existent with the method 250.

Although the marker point identification method 60, the marker pointsub-pixel localization method 110, the physical marker 3D localizationmethod 160, and the marker point mapping method 250 have been describedin connection with CB CT imaging, the methods are also useful for CT,MRI, ultrasound, other forms of medical imaging, and forms ofnon-medical imaging, such as photography. For example, the marker pointidentification method 60 may be used in situations where an object(e.g., a small object) of a well-defined shape needs to be quicklyisolated from a complex background. In addition, although the markerpoint identification method 60 has been described for use with circularor spherically-shaped markers, the method 60 may be used for markerpoints of other well-defined shapes. For example, steps 86 and 88 ofstep 80 (part of method 60) may be modified to recognize other shapes.The method 60 may also use well-defined anatomic objects internal to thepatient as markers. Furthermore, the method 60 may be used as a startingpoint for extracting anatomical landmarks by searching for the generalshape of a particular landmark.

Similarly, the marker point sub-pixel localization method 110 may beapplied to any arbitrarily-shaped marker as long as the characteristicsof the shape are known and the general orientation of the marker in aparticular image can be ascertained. Various approaches can be used toascertain the general orientation of a marker. For example, the markerpoint can be evaluated in general terms using its extensions in one (ortwo) dimensions to help define the placement of the marker point profilein the other dimension. The marker point profile can then be processedwith a simple first-moment calculation or a more elaborate curvefitting. In addition, an image shape can be compared to an atlas ofshapes expected at different orientations and the appropriate one can beselected. Furthermore, reconstruction of CT image data can be used toproduce a 3D image of a marker to determine its orientation beforedefining and fitting a proper curve function.

The marker point sub-pixel localization method 110 also has otherpotential applications in medical imaging and non-medical imaging. Forexample, in dental imaging, the method 110 may be used to place bracketsin an orthodontic treatment, where the positions of the brackets couldbe evaluated to a precision of approximately 100 micrometers or better.The method 110 may also be used for accurate registration of imagesacquired at different times but with the same external or internalmarkers. Similarly, the basic approach of the marker 3D localizationmethod 160 may be used in marker-less CT scanning and other non-CTscanning application, such as astronomy.

The physical marker mapping method 250 also has other potentialapplications than just patient motion tracking. For example, the method250 can be used to create a dynamic calibration for gantries thatexperience non-ideal rotation. In this case, the method 250 would beapplied to a fixed phantom containing embedded markers. This applicationwould be particularly useful for CT systems with x-ray imageintensifiers (sometimes referred to as C-Arm scanning settings).Furthermore, the method 250 can generally be used with systems where anobject has an “ideal” progression that might, in actual practice, not beideal.

In addition, the methods described above can use various data processingtechniques. For example, the marker point mapping method 250 may usevarious lists and combinations of lists or tables for tracking candidatepoints, good points, and suspect points. For example, the method 250 maycreate a candidate point list and remove non-good points rather thancreate a separate good point list. In addition, the method 250 maycreate one or more lists of points and use flags or bits to distinguishbetween candidate points, good points, bad points, and suspect points.

Also, although FIG. 2 illustrates the ROM module 46 as storing separatemodules (e.g., 50, 52, 54, 56, etc.), in some embodiments, these modulesare combined and distributed into one or multiple modules stored in theROM module 46, other memory modules, or combinations thereof Forexample, in some embodiments, the same module performs the markeridentification method and the marker sub-pixel localization method. Inaddition, the ROM module 46 can include other modules for performingfunctions other than those functions described above. In someembodiments, the modules may be stored on a hard-disk, flash, or othersemi-permanent medium and transferred to RAM for execution.

Thus, the invention provides, among other things, methods and systemsfor identifying marker points in an image, determining a sub-pixelcenter point of a marker point identified in an image, determining 3Dcoordinates for a marker identified in an image, and mapping markerpoints to particular physical markers placed on a patient. Variousfeatures and advantages of the invention are set forth in the followingclaims.

1. A method for identifying marker associated points in an image, themethod executed by an imaging system including a scanner, a computerwith an electronic processing unit, and a memory module storing a markerpoint identification module executable by the electronic processingunit, the method comprising: obtaining, at the computer, an image basedon image data generated by the scanner; applying, with the electronicprocessing unit, a two-dimensional high-pass filter with x-raycompensation to the image to suppress background and, thereby, create abackground suppressed image; and extracting, with the electronicprocessing unit, marker associated points from the background suppressedimage.
 2. The method of claim 1, further comprising applying atwo-dimensional low-pass filter to the image to create a low-pass image,and wherein applying the two-dimensional high-pass filter includescombining the image with the low-pass image.
 3. The method of claim 1,wherein applying the two-dimensional high-pass filter includessequentially applying two or more one-dimensional high-pass filters tothe image, each performed in a different two-dimensional direction onthe image.
 4. The method of claim 3, wherein applying eachone-dimensional high-pass filter includes combining at least one of aone-dimensional low-pass filtered representation of the image with theimage, and a one-dimensional low-pass filtered representation of atemporary image from an earlier application of the one-dimensionalhigh-pass filter with the temporary image.
 5. The method of claim 4,wherein the one-dimensional low-pass filtered representation of theimage is created by applying cumulative mean filtering, the filteringcomprising: creating, using the electronic processing unit, a cumulativeimage from the image, wherein each pixel of the cumulative imagerepresents a cumulative sum of lower pixels in a direction of thefilter; shifting, with the electronic processing unit, the cumulativeimage by a filtering factor to create a shifted image; subtracting, withthe electronic processing unit, the cumulative image from the shiftedimage to create a difference image; dividing, with the electronicprocessing unit, the difference image by the filtering factor.
 6. Themethod of claim 3, wherein combining the image and the one-dimensionallow-pass filtered representation of the image accounts for non-linearattenuation of x-rays and includes a pixel-by-pixel division of theimage by the one-dimensional low-pass filtered representation of theimage.
 7. The method of claim 6, wherein the pixel-by-pixel division isperformed using pixels of the one-dimensional low-pass filteredrepresentation of the image as a numerator.
 8. The method of claim 3,wherein combining the image and the one-dimensional low-pass filteredrepresentation of the image includes subtraction.
 9. The method of claim1, wherein extracting the marker associated points from the backgroundsuppressed image includes defining local maxima in the backgroundsuppressed image.
 10. The method of claim 9, wherein defining the localmaxima in the background suppressed image includes identifying pixels inthe background suppressed image that have eight neighboring pixels withsmaller pixel values.
 11. The method of claim 9, wherein extracting themarker associated points from the background suppressed image furtherincludes growing candidate regions based on the local maxima.
 12. Themethod of claim 11, wherein extracting the marker associated points fromthe background suppressed image further includes applying shapecriterion to each candidate region.
 13. The method of claim 12, whereinapplying the shape criterion to each candidate region includesdetermining a ratio of a pixel area of each candidate region and an areadefined by a mean radius of each candidate region squared times pi. 14.The method of claim 11, wherein extracting the marker associated pointsfrom the background suppressed image further includes eliminatingcandidate regions by comparing candidate regions in adjacent projectionimages and applying nearness criteria between candidate regions.
 15. Animaging system including a scanner, a computer with an electronicprocessing unit, and a memory module storing a marker pointidentification module executable by the electronic processing unit,wherein the marker point identification module, when executed, isconfigured to: obtain, at the computer, an image based on image datagenerated by the scanner; apply a two-dimensional high-pass filter withx-ray compensation to the image to suppress background and, thereby,create a background suppressed image; and extract marker associatedpoints from the background suppressed image.
 16. The imaging system ofclaim 15, wherein the marker point identification module, when executed,is further configured to apply a two-dimensional low-pass filter to theimage to create a low-pass image, and apply the two-dimensionalhigh-pass filter by combining the image with the low-pass image.
 17. Theimaging system of claim 15, wherein the marker point identificationmodule, when executed, is further configured to apply thetwo-dimensional high-pass filter by sequentially applying two or moreone-dimensional high-pass filters to the image, each performed in adifferent two-dimensional direction on the image.
 18. The imaging systemof claim 15, wherein the marker point identification module, whenexecuted, is further configured to extract the marker associated pointsfrom the background suppressed image by defining local maxima in thebackground suppressed image, wherein defining local maxima includesidentifying pixels in the background suppressed image that have eightneighboring pixels with smaller pixel values.
 19. The imaging system ofclaim 18, wherein the marker point identification module, when executed,is further configured to extract the marker associated points from thebackground suppressed image by growing candidate regions based on thelocal maxima, applying shape criterion to each candidate region, andeliminating candidate regions by comparing candidate regions in adjacentprojection images and applying nearness criteria between candidateregions.